Data analysis for:

The Foraging-mode Paradigm: A Historical Overview Including a Reevaluation of its Predictions

Dylan Padilla, School of Life Sciences, Arizona State University, Tempe, AZ 85287, United States.

Library

library(AICcmodavg)
library(emo)
library(icons)
library(magick)
library(knitr)
library(MuMIn)
library(png)
library(raster)
library(scales)
library(vioplot)
library(xtable)


Loading data

data <- read.csv("data.csv")
head(data)
>   initial_mass final_mass strain env adult comments
> 1    0.0007245         NA      s   p    NA     died
> 2    0.0007720         NA      s   p    NA     died
> 3    0.0007095         NA      s   u    NA     died
> 4    0.0005566         NA      s   u    NA     died
> 5    0.0008130         NA      s   u    NA     died
> 6    0.0004055  0.0008495      s   p    NA     <NA>
data$growth <-  data$final_mass - data$initial_mass
names(data)
> [1] "initial_mass" "final_mass"   "strain"       "env"          "adult"        "comments"     "growth"
data <- data[ , c(3, 4, 7)]
data <- data[!is.na(data$growth), ]
data <- data[data$growth > 0, ]

table(data$strain, data$env)
>    
>      p  u
>   r 29 29
>   s 32 22
data$strain[data$strain == "s"] <- "Sitter"
data$strain[data$strain == "r"] <- "Rover"
data$env[data$env == "p"] <- "Patchy"
data$env[data$env == "u"] <- "Lumpy"

data <- data[data$growth < 0.002, ] # Removing abnormal data


## Movement data

vid.dat <- read.csv('./OneDrive-2024-08-19/Results/Results_by_ind.csv', row.names = 1, sep = ';')
str(vid.dat)
> 'data.frame': 92 obs. of  25 variables:
>  $ Arena                                : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Sequence                             : chr  "General" "General" "General" "General" ...
>  $ Individual                           : chr  "Ind0" "Ind0" "Ind0" "Ind0" ...
>  $ Beginning_seq                        : num  0 0 0 0 0 0 0 0 0 0 ...
>  $ End_seq                              : num  450 450 450 450 450 ...
>  $ Prop_time_lost                       : num  0 0 0 0 0 0 0 0 0 0 ...
>  $ Smoothing_filter_window              : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Smoothing_Polyorder                  : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Moving_threshold                     : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Prop_time_moving                     : num  0.976 0.949 0.919 0.993 0.945 ...
>  $ Average_Speed                        : num  0.518 2.36 5.962 5.412 2.979 ...
>  $ Average_Speed_Moving                 : num  0.531 2.488 6.486 5.448 3.152 ...
>  $ Traveled_Dist                        : num  233 1062 2683 2435 1341 ...
>  $ Meander                              : num  1423 412 1387 1112 2120 ...
>  $ Traveled_Dist_Moving                 : num  233 1062 2683 2435 1341 ...
>  $ Meander_moving                       : num  1423 412 1387 1112 2120 ...
>  $ Exploration_absolute_value           : num  65.2 1148.3 150.6 1059.7 504.2 ...
>  $ Exploration_relative_value           : num  0.00671 0.12435 0.01568 0.11397 0.04999 ...
>  $ Exploration_method                   : chr  "Modern" "Modern" "Modern" "Modern" ...
>  $ Exploration_area                     : int  10 10 10 10 10 10 10 10 10 10 ...
>  $ Exploration_aspect_param             : logi  NA NA NA NA NA NA ...
>  $ Mean_nb_neighbours                   : logi  NA NA NA NA NA NA ...
>  $ Prop_time_with_at_least_one_neighbour: logi  NA NA NA NA NA NA ...
>  $ Mean_shortest_interind_distance      : logi  NA NA NA NA NA NA ...
>  $ Mean_sum_interind_distances          : logi  NA NA NA NA NA NA ...
vid.dat <- vid.dat[c(2, 4, 1, 3, 5:92), ]
vid.dat[1:8, 1:5]
>                                         Arena Sequence Individual Beginning_seq End_seq
> Video SP - Apr 14 2024, 12 14 18.avi(3)     0  General       Ind0             0  449.97
> Video RP - Apr 14 2024, 12 14 18.avi(4)     0  General       Ind0             0  449.97
> Video SU - Apr 14 2024, 12 14 18.avi(1)     0  General       Ind0             0  449.97
> Video RU - Apr 14 2024, 12 14 18.avi(2)     0  General       Ind0             0  449.97
> Video SP - Apr 14 2024, 14 03 52.avi(1)     0  General       Ind0             0  449.97
> Video RP - Apr 14 2024, 14 03 52.avi(2)     0  General       Ind0             0  449.97
> Video SU - Apr 14 2024, 14 03 52.avi(3)     0  General       Ind0             0  449.97
> Video RU - Apr 14 2024, 14 03 52.avi(4)     0  General       Ind0             0  449.97
mov.dat <- read.csv('movement.csv')
str(mov.dat)
> 'data.frame': 92 obs. of  6 variables:
>  $ mass    : num  0.001991 0.001548 0.000821 0.001388 0.001131 ...
>  $ strain  : chr  "r" "r" "s" "s" ...
>  $ env     : chr  "p" "u" "p" "u" ...
>  $ dis_trav: logi  NA NA NA NA NA NA ...
>  $ video   : int  1 1 1 1 2 2 2 2 3 3 ...
>  $ comments: logi  NA NA NA NA NA NA ...
head(mov.dat)
>        mass strain env dis_trav video comments
> 1 0.0019910      r   p       NA     1       NA
> 2 0.0015480      r   u       NA     1       NA
> 3 0.0008210      s   p       NA     1       NA
> 4 0.0013880      s   u       NA     1       NA
> 5 0.0011315      r   p       NA     2       NA
> 6 0.0014595      r   u       NA     2       NA
s <- mov.dat[mov.dat$strain == "s", ]
r <- mov.dat[mov.dat$strain == "r", ]
dim(s)
> [1] 46  6
dim(r)
> [1] 46  6
ord.dat <- data.frame()

for(i in 1:46){
    idx <- s[i, ]
    idx2 <- r[i, ]
    obj <- rbind(idx, idx2)
    ord.dat <- rbind(ord.dat, obj)
}

ord.dat[1:8, ]
>        mass strain env dis_trav video comments
> 3 0.0008210      s   p       NA     1       NA
> 1 0.0019910      r   p       NA     1       NA
> 4 0.0013880      s   u       NA     1       NA
> 2 0.0015480      r   u       NA     1       NA
> 7 0.0014475      s   p       NA     2       NA
> 5 0.0011315      r   p       NA     2       NA
> 8 0.0012095      s   u       NA     2       NA
> 6 0.0014595      r   u       NA     2       NA
vid.dat[1:8, 1:5]
>                                         Arena Sequence Individual Beginning_seq End_seq
> Video SP - Apr 14 2024, 12 14 18.avi(3)     0  General       Ind0             0  449.97
> Video RP - Apr 14 2024, 12 14 18.avi(4)     0  General       Ind0             0  449.97
> Video SU - Apr 14 2024, 12 14 18.avi(1)     0  General       Ind0             0  449.97
> Video RU - Apr 14 2024, 12 14 18.avi(2)     0  General       Ind0             0  449.97
> Video SP - Apr 14 2024, 14 03 52.avi(1)     0  General       Ind0             0  449.97
> Video RP - Apr 14 2024, 14 03 52.avi(2)     0  General       Ind0             0  449.97
> Video SU - Apr 14 2024, 14 03 52.avi(3)     0  General       Ind0             0  449.97
> Video RU - Apr 14 2024, 14 03 52.avi(4)     0  General       Ind0             0  449.97
dim(ord.dat)
> [1] 92  6
vid.dat$mass <- ord.dat$mass
vid.dat$strain <- ord.dat$strain
vid.dat$env <- ord.dat$env

## Proportion of area visited

area.dat <- read.csv('Proportion-of-area-visited.csv', sep = ';')
area.dat <- area.dat[c(2, 4, 1, 3, 5:92), ]
vid.dat$prop_area_visited <- area.dat[ , 2]

names(vid.dat)
>  [1] "Arena"                                 "Sequence"                              "Individual"                           
>  [4] "Beginning_seq"                         "End_seq"                               "Prop_time_lost"                       
>  [7] "Smoothing_filter_window"               "Smoothing_Polyorder"                   "Moving_threshold"                     
> [10] "Prop_time_moving"                      "Average_Speed"                         "Average_Speed_Moving"                 
> [13] "Traveled_Dist"                         "Meander"                               "Traveled_Dist_Moving"                 
> [16] "Meander_moving"                        "Exploration_absolute_value"            "Exploration_relative_value"           
> [19] "Exploration_method"                    "Exploration_area"                      "Exploration_aspect_param"             
> [22] "Mean_nb_neighbours"                    "Prop_time_with_at_least_one_neighbour" "Mean_shortest_interind_distance"      
> [25] "Mean_sum_interind_distances"           "mass"                                  "strain"                               
> [28] "env"                                   "prop_area_visited"
str(vid.dat)
> 'data.frame': 92 obs. of  29 variables:
>  $ Arena                                : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Sequence                             : chr  "General" "General" "General" "General" ...
>  $ Individual                           : chr  "Ind0" "Ind0" "Ind0" "Ind0" ...
>  $ Beginning_seq                        : num  0 0 0 0 0 0 0 0 0 0 ...
>  $ End_seq                              : num  450 450 450 450 450 ...
>  $ Prop_time_lost                       : num  0 0 0 0 0 0 0 0 0 0 ...
>  $ Smoothing_filter_window              : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Smoothing_Polyorder                  : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Moving_threshold                     : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Prop_time_moving                     : num  0.949 0.993 0.976 0.919 0.945 ...
>  $ Average_Speed                        : num  2.36 5.412 0.518 5.962 2.979 ...
>  $ Average_Speed_Moving                 : num  2.488 5.448 0.531 6.486 3.152 ...
>  $ Traveled_Dist                        : num  1062 2435 233 2683 1341 ...
>  $ Meander                              : num  412 1112 1423 1387 2120 ...
>  $ Traveled_Dist_Moving                 : num  1062 2435 233 2683 1341 ...
>  $ Meander_moving                       : num  412 1112 1423 1387 2120 ...
>  $ Exploration_absolute_value           : num  1148.3 1059.7 65.2 150.6 504.2 ...
>  $ Exploration_relative_value           : num  0.12435 0.11397 0.00671 0.01568 0.04999 ...
>  $ Exploration_method                   : chr  "Modern" "Modern" "Modern" "Modern" ...
>  $ Exploration_area                     : int  10 10 10 10 10 10 10 10 10 10 ...
>  $ Exploration_aspect_param             : logi  NA NA NA NA NA NA ...
>  $ Mean_nb_neighbours                   : logi  NA NA NA NA NA NA ...
>  $ Prop_time_with_at_least_one_neighbour: logi  NA NA NA NA NA NA ...
>  $ Mean_shortest_interind_distance      : logi  NA NA NA NA NA NA ...
>  $ Mean_sum_interind_distances          : logi  NA NA NA NA NA NA ...
>  $ mass                                 : num  0.000821 0.001991 0.001388 0.001548 0.001448 ...
>  $ strain                               : chr  "s" "r" "s" "r" ...
>  $ env                                  : chr  "p" "p" "u" "u" ...
>  $ prop_area_visited                    : num  0.124 0.114 0.007 0.016 0.05 0.147 0.003 0.041 0.24 0.17 ...
data.frame(vid.dat[1:12, 10:11])
>                                         Prop_time_moving Average_Speed
> Video SP - Apr 14 2024, 12 14 18.avi(3)        0.9486629     2.3603480
> Video RP - Apr 14 2024, 12 14 18.avi(4)        0.9934810     5.4120290
> Video SU - Apr 14 2024, 12 14 18.avi(1)        0.9757760     0.5177664
> Video RU - Apr 14 2024, 12 14 18.avi(2)        0.9191051     5.9617525
> Video SP - Apr 14 2024, 14 03 52.avi(1)        0.9452552     2.9791838
> Video RP - Apr 14 2024, 14 03 52.avi(2)        0.9984443     1.2822174
> Video SU - Apr 14 2024, 14 03 52.avi(3)        0.9691088     0.3256189
> Video RU - Apr 14 2024, 14 03 52.avi(4)        0.9758501     5.5721335
> Video SP - Apr 14 2024, 15 13 46.avi(1)        0.9748870     2.7576846
> Video RP - Apr 14 2024, 15 13 46.avi(2)        0.9659234     1.3956371
> Video SU - Apr 14 2024, 15 13 46.avi(3)        0.9698496     2.8616386
> Video RU - Apr 14 2024, 15 13 46.avi(4)        0.9392548     3.1724128


Modeling foraging behavior and growth

##vid.dat <- vid.dat[vid.dat$Prop_time_moving > 0.7, ] ## extreme value
vid.dat$env[vid.dat$env == 'u'] <- 'Lumpy'
vid.dat$env[vid.dat$env == 'p'] <- 'Patchy'
vid.dat$strain[vid.dat$strain == 'r'] <- 'Rover'
vid.dat$strain[vid.dat$strain == 's'] <- 'Sitter'
str(vid.dat)
> 'data.frame': 92 obs. of  29 variables:
>  $ Arena                                : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Sequence                             : chr  "General" "General" "General" "General" ...
>  $ Individual                           : chr  "Ind0" "Ind0" "Ind0" "Ind0" ...
>  $ Beginning_seq                        : num  0 0 0 0 0 0 0 0 0 0 ...
>  $ End_seq                              : num  450 450 450 450 450 ...
>  $ Prop_time_lost                       : num  0 0 0 0 0 0 0 0 0 0 ...
>  $ Smoothing_filter_window              : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Smoothing_Polyorder                  : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Moving_threshold                     : int  0 0 0 0 0 0 0 0 0 0 ...
>  $ Prop_time_moving                     : num  0.949 0.993 0.976 0.919 0.945 ...
>  $ Average_Speed                        : num  2.36 5.412 0.518 5.962 2.979 ...
>  $ Average_Speed_Moving                 : num  2.488 5.448 0.531 6.486 3.152 ...
>  $ Traveled_Dist                        : num  1062 2435 233 2683 1341 ...
>  $ Meander                              : num  412 1112 1423 1387 2120 ...
>  $ Traveled_Dist_Moving                 : num  1062 2435 233 2683 1341 ...
>  $ Meander_moving                       : num  412 1112 1423 1387 2120 ...
>  $ Exploration_absolute_value           : num  1148.3 1059.7 65.2 150.6 504.2 ...
>  $ Exploration_relative_value           : num  0.12435 0.11397 0.00671 0.01568 0.04999 ...
>  $ Exploration_method                   : chr  "Modern" "Modern" "Modern" "Modern" ...
>  $ Exploration_area                     : int  10 10 10 10 10 10 10 10 10 10 ...
>  $ Exploration_aspect_param             : logi  NA NA NA NA NA NA ...
>  $ Mean_nb_neighbours                   : logi  NA NA NA NA NA NA ...
>  $ Prop_time_with_at_least_one_neighbour: logi  NA NA NA NA NA NA ...
>  $ Mean_shortest_interind_distance      : logi  NA NA NA NA NA NA ...
>  $ Mean_sum_interind_distances          : logi  NA NA NA NA NA NA ...
>  $ mass                                 : num  0.000821 0.001991 0.001388 0.001548 0.001448 ...
>  $ strain                               : chr  "Sitter" "Rover" "Sitter" "Rover" ...
>  $ env                                  : chr  "Patchy" "Patchy" "Lumpy" "Lumpy" ...
>  $ prop_area_visited                    : num  0.124 0.114 0.007 0.016 0.05 0.147 0.003 0.041 0.24 0.17 ...
## Proportion of area covered


mod <- lm(prop_area_visited ~ mass*strain*env, data = vid.dat)
summary(mod)
> 
> Call:
> lm(formula = prop_area_visited ~ mass * strain * env, data = vid.dat)
> 
> Residuals:
>      Min       1Q   Median       3Q      Max 
> -0.15443 -0.06215 -0.02083  0.06943  0.21500 
> 
> Coefficients:
>                               Estimate Std. Error t value Pr(>|t|)  
> (Intercept)                    0.11723    0.07979   1.469   0.1455  
> mass                          18.15868   54.18088   0.335   0.7383  
> strainSitter                  -0.20760    0.12371  -1.678   0.0970 .
> envPatchy                      0.28803    0.13347   2.158   0.0338 *
> mass:strainSitter             95.17618   88.27340   1.078   0.2840  
> mass:envPatchy              -152.05488   95.00854  -1.600   0.1133  
> strainSitter:envPatchy        -0.04573    0.18349  -0.249   0.8038  
> mass:strainSitter:envPatchy    6.03000  134.26488   0.045   0.9643  
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Residual standard error: 0.09032 on 84 degrees of freedom
> Multiple R-squared:  0.3546,  Adjusted R-squared:  0.3008 
> F-statistic: 6.592 on 7 and 84 DF,  p-value: 3.273e-06
mod1 <- lm(prop_area_visited ~ mass*strain + env, data = vid.dat)
summary(mod1)
> 
> Call:
> lm(formula = prop_area_visited ~ mass * strain + env, data = vid.dat)
> 
> Residuals:
>      Min       1Q   Median       3Q      Max 
> -0.15481 -0.06221 -0.02851  0.07036  0.21089 
> 
> Coefficients:
>                    Estimate Std. Error t value Pr(>|t|)    
> (Intercept)         0.19850    0.06548   3.031 0.003209 ** 
> mass              -34.18583   44.90920  -0.761 0.448584    
> strainSitter       -0.19425    0.08952  -2.170 0.032746 *  
> envPatchy           0.06584    0.01925   3.420 0.000956 ***
> mass:strainSitter  71.56717   65.26893   1.096 0.275887    
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Residual standard error: 0.09156 on 87 degrees of freedom
> Multiple R-squared:  0.3131,  Adjusted R-squared:  0.2815 
> F-statistic: 9.915 on 4 and 87 DF,  p-value: 1.173e-06
mod2 <- lm(prop_area_visited ~ strain*env, data = vid.dat)
summary(mod2)
> 
> Call:
> lm(formula = prop_area_visited ~ strain * env, data = vid.dat)
> 
> Residuals:
>      Min       1Q   Median       3Q      Max 
> -0.13956 -0.06178 -0.03193  0.07027  0.21578 
> 
> Coefficients:
>                        Estimate Std. Error t value Pr(>|t|)    
> (Intercept)             0.14322    0.01904   7.521 4.36e-11 ***
> strainSitter           -0.08291    0.02693  -3.079  0.00277 ** 
> envPatchy               0.08135    0.02693   3.021  0.00330 ** 
> strainSitter:envPatchy -0.03087    0.03808  -0.811  0.41981    
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Residual standard error: 0.09132 on 88 degrees of freedom
> Multiple R-squared:  0.3088,  Adjusted R-squared:  0.2852 
> F-statistic:  13.1 on 3 and 88 DF,  p-value: 3.767e-07
anova(mod2)
> Analysis of Variance Table
> 
> Response: prop_area_visited
>            Df  Sum Sq  Mean Sq F value    Pr(>F)    
> strain      1 0.22246 0.222463  26.674 1.483e-06 ***
> env         1 0.09992 0.099924  11.981 0.0008313 ***
> strain:env  1 0.00548 0.005479   0.657 0.4198067    
> Residuals  88 0.73391 0.008340                      
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod3 <- lm(prop_area_visited ~ strain + env, data = vid.dat)
summary(mod3)
> 
> Call:
> lm(formula = prop_area_visited ~ strain + env, data = vid.dat)
> 
> Residuals:
>      Min       1Q   Median       3Q      Max 
> -0.14694 -0.06400 -0.02426  0.07215  0.20806 
> 
> Coefficients:
>              Estimate Std. Error t value Pr(>|t|)    
> (Intercept)   0.15093    0.01646   9.170 1.67e-14 ***
> strainSitter -0.09835    0.01901  -5.175 1.40e-06 ***
> envPatchy     0.06591    0.01901   3.468  0.00081 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Residual standard error: 0.09115 on 89 degrees of freedom
> Multiple R-squared:  0.3036,  Adjusted R-squared:  0.288 
> F-statistic:  19.4 on 2 and 89 DF,  p-value: 1.015e-07
anova(mod3)
> Analysis of Variance Table
> 
> Response: prop_area_visited
>           Df  Sum Sq  Mean Sq F value    Pr(>F)    
> strain     1 0.22246 0.222463  26.778 1.399e-06 ***
> env        1 0.09992 0.099924  12.028 0.0008096 ***
> Residuals 89 0.73939 0.008308                      
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Model selection procedure based on AIC criterion

output.models <- model.sel(mod1, mod2, mod3)

sel.table <- round(as.data.frame(output.models)[7:11], 3)
names(sel.table)[1] <- "K"
sel.table$Model <- rownames(sel.table)
sel.table <- sel.table[ , c(6, 1, 2, 3, 4, 5)]
sel.table
>      Model K logLik     AICc delta weight
> mod3  mod3 4 91.349 -174.237 0.000  0.604
> mod2  mod2 5 91.691 -172.684 1.554  0.278
> mod1  mod1 6 91.980 -170.972 3.265  0.118
xtable(sel.table, digits = 3, caption = 'Caption here')
> % latex table generated in R 4.4.1 by xtable 1.8-4 package
> % Sat Aug 24 10:49:20 2024
> \begin{table}[ht]
> \centering
> \begin{tabular}{rlrrrrr}
>   \hline
>  & Model & K & logLik & AICc & delta & weight \\ 
>   \hline
> mod3 & mod3 & 4.000 & 91.349 & -174.237 & 0.000 & 0.604 \\ 
>   mod2 & mod2 & 5.000 & 91.691 & -172.684 & 1.554 & 0.278 \\ 
>   mod1 & mod1 & 6.000 & 91.980 & -170.972 & 3.265 & 0.118 \\ 
>    \hline
> \end{tabular}
> \caption{Caption here} 
> \end{table}
## Model diagnosis

layout(matrix(c(0, 0, 0, 0,
                1, 1, 2, 2,
                1, 1, 2, 2,
                0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))



## Checking homogeneity of variance


plot(fitted(mod3), resid(mod3), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", type = 'n', las = 1)
grid()

par(new = TRUE)
plot(fitted(mod3), resid(mod3), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", las = 1)

abline(h = 0, col = "darkorange", lwd = 2)


## Checking normality

qqnorm(resid(mod3), col = "darkgrey", type = 'n', las = 1)
grid()

par(new = TRUE)
qqnorm(resid(mod3), col = "darkgrey", las = 1)
qqline(resid(mod3), col = "dodgerblue", lwd = 2)

## Distance traveled

mod4 <- lm(Traveled_Dist_Moving ~ mass*strain + env, data = vid.dat)
summary(mod4)
> 
> Call:
> lm(formula = Traveled_Dist_Moving ~ mass * strain + env, data = vid.dat)
> 
> Residuals:
>      Min       1Q   Median       3Q      Max 
> -1009.14  -424.53   -53.04   350.55  1555.33 
> 
> Coefficients:
>                    Estimate Std. Error t value Pr(>|t|)  
> (Intercept)           566.2      396.1   1.429   0.1565  
> mass               362471.5   271648.8   1.334   0.1856  
> strainSitter         -169.6      541.5  -0.313   0.7549  
> envPatchy             210.1      116.5   1.804   0.0747 .
> mass:strainSitter -287265.1   394801.6  -0.728   0.4688  
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Residual standard error: 553.8 on 87 degrees of freedom
> Multiple R-squared:  0.2534,  Adjusted R-squared:  0.2191 
> F-statistic: 7.381 on 4 and 87 DF,  p-value: 3.629e-05
mod5 <- lm(Traveled_Dist_Moving ~ strain*env, data = vid.dat)
summary(mod5)
> 
> Call:
> lm(formula = Traveled_Dist_Moving ~ strain * env, data = vid.dat)
> 
> Residuals:
>     Min      1Q  Median      3Q     Max 
> -1130.4  -411.9  -105.2   299.3  1534.0 
> 
> Coefficients:
>                        Estimate Std. Error t value Pr(>|t|)    
> (Intercept)             1148.55     115.07   9.981 3.93e-16 ***
> strainSitter            -715.72     162.74  -4.398 3.05e-05 ***
> envPatchy                 53.19     162.74   0.327    0.745    
> strainSitter:envPatchy   278.99     230.15   1.212    0.229    
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Residual standard error: 551.9 on 88 degrees of freedom
> Multiple R-squared:  0.2501,  Adjusted R-squared:  0.2245 
> F-statistic: 9.783 on 3 and 88 DF,  p-value: 1.233e-05
anova(mod5)
> Analysis of Variance Table
> 
> Response: Traveled_Dist_Moving
>            Df   Sum Sq Mean Sq F value    Pr(>F)    
> strain      1  7636784 7636784 25.0742 2.813e-06 ***
> env         1   853966  853966  2.8039   0.09759 .  
> strain:env  1   447556  447556  1.4695   0.22867    
> Residuals  88 26801948  304568                      
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod6 <- lm(Traveled_Dist_Moving ~ strain + env, data = vid.dat)
summary(mod6)
> 
> Call:
> lm(formula = Traveled_Dist_Moving ~ strain + env, data = vid.dat)
> 
> Residuals:
>     Min      1Q  Median      3Q     Max 
> -1060.7  -414.0  -109.2   366.2  1603.8 
> 
> Coefficients:
>              Estimate Std. Error t value Pr(>|t|)    
> (Intercept)   1078.80      99.92  10.797  < 2e-16 ***
> strainSitter  -576.22     115.38  -4.994 2.92e-06 ***
> envPatchy      192.69     115.38   1.670   0.0984 .  
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Residual standard error: 553.3 on 89 degrees of freedom
> Multiple R-squared:  0.2376,  Adjusted R-squared:  0.2204 
> F-statistic: 13.87 on 2 and 89 DF,  p-value: 5.727e-06
anova(mod6)
> Analysis of Variance Table
> 
> Response: Traveled_Dist_Moving
>           Df   Sum Sq Mean Sq F value    Pr(>F)    
> strain     1  7636784 7636784 24.9426 2.921e-06 ***
> env        1   853966  853966  2.7891   0.09842 .  
> Residuals 89 27249504  306174                      
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod7 <- lm(Traveled_Dist_Moving ~ strain, data = vid.dat)
summary(mod7)
> 
> Call:
> lm(formula = Traveled_Dist_Moving ~ strain, data = vid.dat)
> 
> Residuals:
>      Min       1Q   Median       3Q      Max 
> -1157.02  -495.06   -59.71   351.98  1507.44 
> 
> Coefficients:
>              Estimate Std. Error t value Pr(>|t|)    
> (Intercept)   1175.15      82.39  14.263  < 2e-16 ***
> strainSitter  -576.22     116.52  -4.945 3.51e-06 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Residual standard error: 558.8 on 90 degrees of freedom
> Multiple R-squared:  0.2137,  Adjusted R-squared:  0.2049 
> F-statistic: 24.46 on 1 and 90 DF,  p-value: 3.507e-06
anova(mod7)
> Analysis of Variance Table
> 
> Response: Traveled_Dist_Moving
>           Df   Sum Sq Mean Sq F value    Pr(>F)    
> strain     1  7636784 7636784  24.456 3.507e-06 ***
> Residuals 90 28103469  312261                      
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Model selection procedure based on AIC criterion

output.models <- model.sel(mod4, mod5, mod6, mod7)

sel.table <- round(as.data.frame(output.models)[7:11], 3)
names(sel.table)[1] <- "K"
sel.table$Model <- rownames(sel.table)
sel.table <- sel.table[ , c(6, 1, 2, 3, 4, 5)]
sel.table
>      Model K   logLik     AICc delta weight
> mod6  mod6 4 -710.085 1428.630 0.000  0.371
> mod7  mod7 3 -711.505 1429.282 0.652  0.268
> mod5  mod5 5 -709.323 1429.344 0.714  0.260
> mod4  mod4 6 -709.121 1431.231 2.600  0.101
## Model diagnosis


layout(matrix(c(0, 0, 0, 0,
                1, 1, 2, 2,
                1, 1, 2, 2,
                0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))



## Checking homogeneity of variance


plot(fitted(mod6), resid(mod6), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", type = 'n', las = 1)
grid()

par(new = TRUE)
plot(fitted(mod6), resid(mod6), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", las = 1)

abline(h = 0, col = "darkorange", lwd = 2)


## Checking normality

qqnorm(resid(mod6), col = "darkgrey", type = 'n', las = 1)
grid()

par(new = TRUE)
qqnorm(resid(mod6), col = "darkgrey", las = 1)
qqline(resid(mod6), col = "dodgerblue", lwd = 2)

## Growth rate

int.model <- lm(growth ~ strain*env, data = data)
summary(int.model)
> 
> Call:
> lm(formula = growth ~ strain * env, data = data)
> 
> Residuals:
>        Min         1Q     Median         3Q        Max 
> -3.238e-04 -1.061e-04  1.380e-06  1.033e-04  5.414e-04 
> 
> Coefficients:
>                          Estimate Std. Error t value Pr(>|t|)    
> (Intercept)             3.221e-04  3.001e-05  10.735  < 2e-16 ***
> strainSitter            1.522e-05  4.569e-05   0.333 0.739690    
> envPatchy              -1.655e-04  4.321e-05  -3.830 0.000218 ***
> strainSitter:envPatchy  3.126e-05  6.242e-05   0.501 0.617629    
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Residual standard error: 0.0001616 on 105 degrees of freedom
> Multiple R-squared:  0.185,   Adjusted R-squared:  0.1618 
> F-statistic: 7.947 on 3 and 105 DF,  p-value: 7.955e-05
anova(int.model)
> Analysis of Variance Table
> 
> Response: growth
>             Df     Sum Sq    Mean Sq F value    Pr(>F)    
> strain       1 7.4100e-09 7.4100e-09  0.2836    0.5955    
> env          1 6.0859e-07 6.0859e-07 23.3073 4.712e-06 ***
> strain:env   1 6.5500e-09 6.5500e-09  0.2507    0.6176    
> Residuals  105 2.7417e-06 2.6110e-08                      
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
noint.model <- lm(growth ~ strain + env, data = data)
summary(noint.model)
> 
> Call:
> lm(formula = growth ~ strain + env, data = data)
> 
> Residuals:
>        Min         1Q     Median         3Q        Max 
> -0.0003334 -0.0001039  0.0000036  0.0001027  0.0005486 
> 
> Coefficients:
>                Estimate Std. Error t value Pr(>|t|)    
> (Intercept)   3.149e-04  2.622e-05  12.010  < 2e-16 ***
> strainSitter  3.196e-05  3.102e-05   1.030    0.305    
> envPatchy    -1.506e-04  3.107e-05  -4.845 4.35e-06 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Residual standard error: 0.000161 on 106 degrees of freedom
> Multiple R-squared:  0.1831,  Adjusted R-squared:  0.1677 
> F-statistic: 11.88 on 2 and 106 DF,  p-value: 2.213e-05
anova(noint.model)
> Analysis of Variance Table
> 
> Response: growth
>            Df     Sum Sq    Mean Sq F value    Pr(>F)    
> strain      1 7.4100e-09 7.4100e-09  0.2856    0.5942    
> env         1 6.0859e-07 6.0859e-07 23.4732 4.347e-06 ***
> Residuals 106 2.7482e-06 2.5930e-08                      
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Model selection procedure based on AIC criterion

output.models <- model.sel(int.model, noint.model)

sel.table <- round(as.data.frame(output.models)[5:9], 3)
names(sel.table)[1] <- "K"
sel.table$Model <- rownames(sel.table)
sel.table <- sel.table[ , c(6, 1, 2, 3, 4, 5)]
sel.table
>                   Model K  logLik      AICc delta weight
> noint.model noint.model 4 798.862 -1589.340 0.000  0.725
> int.model     int.model 5 798.992 -1587.402 1.938  0.275
xtable(sel.table, digits = 3, caption = 'Caption here')
> % latex table generated in R 4.4.1 by xtable 1.8-4 package
> % Sat Aug 24 10:49:20 2024
> \begin{table}[ht]
> \centering
> \begin{tabular}{rlrrrrr}
>   \hline
>  & Model & K & logLik & AICc & delta & weight \\ 
>   \hline
> noint.model & noint.model & 4.000 & 798.862 & -1589.340 & 0.000 & 0.725 \\ 
>   int.model & int.model & 5.000 & 798.992 & -1587.402 & 1.938 & 0.275 \\ 
>    \hline
> \end{tabular}
> \caption{Caption here} 
> \end{table}
## Model diagnosis


layout(matrix(c(0, 0, 0, 0,
                1, 1, 2, 2,
                1, 1, 2, 2,
                0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))



## Checking homogeneity of variance


plot(fitted(noint.model), resid(noint.model), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", type = 'n', las = 1)
grid()

par(new = TRUE)
plot(fitted(noint.model), resid(noint.model), col = "grey", pch = 20, xlab = "Fitted", ylab = "Residual", main = "Fitted versus Residuals", las = 1)

abline(h = 0, col = "darkorange", lwd = 2)


## Checking normality

qqnorm(resid(noint.model), col = "darkgrey", type = 'n', las = 1)
grid()

par(new = TRUE)
qqnorm(resid(noint.model), col = "darkgrey", las = 1)
qqline(resid(noint.model), col = "dodgerblue", lwd = 2)

## Proportion of area covered

xtable(data.frame(anova((mod3))), caption = 'Caption here...', digits = 3)
> % latex table generated in R 4.4.1 by xtable 1.8-4 package
> % Sat Aug 24 10:49:21 2024
> \begin{table}[ht]
> \centering
> \begin{tabular}{rrrrrr}
>   \hline
>  & Df & Sum.Sq & Mean.Sq & F.value & Pr..F. \\ 
>   \hline
> strain &    1 & 0.222 & 0.222 & 26.778 & 0.000 \\ 
>   env &    1 & 0.100 & 0.100 & 12.028 & 0.001 \\ 
>   Residuals &   89 & 0.739 & 0.008 &  &  \\ 
>    \hline
> \end{tabular}
> \caption{Caption here...} 
> \end{table}
## Growth

xtable(formatC(as.matrix(anova(noint.model)), digits = 3, format = 'e'), caption = 'Caption here...', digits = 3)
> % latex table generated in R 4.4.1 by xtable 1.8-4 package
> % Sat Aug 24 10:49:21 2024
> \begin{table}[ht]
> \centering
> \begin{tabular}{rlllll}
>   \hline
>  & Df & Sum Sq & Mean Sq & F value & Pr($>$F) \\ 
>   \hline
> strain & 1.000e+00 & 7.405e-09 & 7.405e-09 & 2.856e-01 & 5.942e-01 \\ 
>   env & 1.000e+00 & 6.086e-07 & 6.086e-07 & 2.347e+01 & 4.347e-06 \\ 
>   Residuals & 1.060e+02 & 2.748e-06 & 2.593e-08 &   NA &   NA \\ 
>    \hline
> \end{tabular}
> \caption{Caption here...} 
> \end{table}


Plotting results of the model

## Barplot

## png('figure2.png', height = 7, width = 7, units = 'in', res = 350)

## Let's calculate the average value for each condition and each specie with the 'tapply' function

mean.area <- tapply(vid.dat$prop_area_visited, list(vid.dat$strain, vid.dat$env), mean)
as.matrix(mean.area)
>             Lumpy    Patchy
> Rover  0.14321739 0.2245652
> Sitter 0.06030435 0.1107826
## Plot boundaries

lim <- 2*max(mean.area)

## A function to add arrows on the chart

error.bar <- function(x, y, upper, lower = upper, length = 0.1,...){
  arrows(x, y+upper, x, y-lower, angle = 90, code = 3, length = length, ...)
}
 
## Then I calculate the standard deviation for each specie and condition :

std.area <- tapply(vid.dat$prop_area_visited, list(vid.dat$strain, vid.dat$env), sd)
std.area
>             Lumpy     Patchy
> Rover  0.10867171 0.07829944
> Sitter 0.07042819 0.10227012
stdev.area <- as.matrix(std.area*1.96/46)

## I am ready to add the error bar on the plot using my "error bar" function!

cols <- c("black" , "gray")

par(mar = c(5, 6, 1, 1), mgp = c(3.8, 1, 0))

ze_barplot <- barplot(NA, beside = TRUE, legend.text = T, col = alpha(cols, 0.2), ylim = c(0, lim), ylab = expression(paste('Proportion of area visited')~(mm^2)), xlab = "Environment", las = 1)
grid()

par(new = TRUE)
ze_barplot <- barplot(mean.area, beside = TRUE, legend.text = T, col = alpha(cols, 0.5) , ylim = c(0, lim), ylab = expression(paste('Proportion of area visited')~(mm^2)), xlab = "Environment", las = 1)
error.bar(ze_barplot, mean.area, stdev.area, col = "black")
box()

## dev.off()
## Barplot

## png('figure2.png', height = 7, width = 7, units = 'in', res = 350)

## Let's calculate the average value for each condition and each specie with the 'tapply' function

mean.gr <- tapply(data$growth, list(data$strain, data$env), mean)
as.matrix(mean.gr)
>               Lumpy       Patchy
> Rover  0.0003221207 0.0001565889
> Sitter 0.0003373409 0.0002030645
## Plot boundaries

lim <- 2*max(mean.gr)

## A function to add arrows on the chart

error.bar <- function(x, y, upper, lower = upper, length = 0.1,...){
  arrows(x, y+upper, x, y-lower, angle = 90, code = 3, length = length, ...)
}
 
## Then I calculate the standard deviation for each specie and condition :

std.gr <- tapply(data$growth, list(data$strain, data$env), sd)
std.gr
>               Lumpy       Patchy
> Rover  0.0001965322 9.225672e-05
> Sitter 0.0001966897 1.445087e-04
stdev <- as.matrix(std.gr*1.96/46)

## I am ready to add the error bar on the plot using my "error bar" function!

cols <- c("black" , "gray")

par(mar = c(5, 6, 1, 1), mgp = c(3.8, 1, 0))

ze_barplot <- barplot(NA, beside = TRUE, legend.text = T, col = alpha(cols, 0.2), ylim = c(0, lim), ylab = expression(Growth~(Delta*gr%*%day^-1)), xlab = "Environment", las = 1)
grid()

par(new = TRUE)
ze_barplot <- barplot(mean.gr, beside = TRUE, legend.text = T, col = alpha(cols, 0.5) , ylim = c(0, lim), ylab = expression(Growth~(Delta*gr%*%day^-1)), xlab = "Environment", las = 1)
error.bar(ze_barplot, mean.gr, stdev, col = "black")
box()

## dev.off()


Figure 1

for(i in 1:23){
    
    png(paste('figure', i, '.png', sep = ''), width = 7, height = 7, units = 'in', res = 350)
    
    layout(matrix(c(1, 1, 2, 2,
                    1, 1, 2, 2,
                    3, 3, 4, 4,
                    3, 3, 4, 4), nrow = 4, ncol = 4, byrow = TRUE))
    
    RP <- readPNG(source = paste('imgs/RP-path/figure', i, '.png', sep = ''))
    SP <- readPNG(source = paste('imgs/SP-path/figure', i, '.png', sep = ''))
    RU <- readPNG(source = paste('imgs/RU-path/figure', i, '.png', sep = ''))
    SU <- readPNG(source = paste('imgs/SU-path/figure', i, '.png', sep = ''))
    play <- readPNG('imgs/play-button.png')
    
    par(mar = c(2.5, 1, 1, 0))
    plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '', main = 'Rovers')
    par(new = TRUE)
    rasterImage(RP, 0, 0, 1, 1)
    mtext('A', at = 0, line = -0.5)

    par(mar = c(2.5, 1, 1, 0))
    plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '', main = 'Sitters')
    par(new = TRUE)
    rasterImage(SP, 0, 0, 1, 1)

    par(mar = c(2.5, 1, 1, 0))
    plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '')
    par(new = TRUE)
    rasterImage(RU, 0, 0, 1, 1)
    mtext('B', at = 0, line = -0.5)

   
    par(mar = c(2.5, 1, 1, 0))
    plot(NA, ylim = c(0, 1), xlim = c(0, 1), type = 'n', axes = FALSE, xlab = '', ylab = '')
    par(new = TRUE)
    rasterImage(SU, 0, 0, 1, 1)

    par(xpd = TRUE)
    rasterImage(play, -0.08, 0, 0.02, 0.11)
    
    dev.off()
}



Figure 2

##png('figure2.png', width = 7, height = 7, units = 'in', res = 350)


layout(matrix(c(0, 1, 1, 0,
                0, 1, 1, 0,
                0, 2, 2, 0,
                0, 2, 2, 0), nrow = 4, ncol = 4, byrow = TRUE))


par(mar = c(5, 6, 1, 1), mgp = c(3.8, 1, 0))

lim <- 2*max(mean.area)

ze_barplot <- barplot(NA, beside = TRUE, legend.text = T, col = alpha(cols, 0.2), ylim = c(0, lim), ylab = expression(paste('Proportion of area visited')~(mm^2)), xlab = "Environment", las = 1)
grid()

par(new = TRUE)
ze_barplot <- barplot(mean.area, beside = TRUE, legend.text = T, col = alpha(cols, 0.5) , ylim = c(0, lim), ylab = expression(paste('Proportion of area visited')~(mm^2)), xlab = "Environment", las = 1)
error.bar(ze_barplot, mean.area, stdev.area, col = "black")
box()
mtext('A', at = -0.45, line = -0.5)


par(mar = c(5, 6, 1, 1), mgp = c(3.8, 1, 0))

lim <- 2*max(mean.gr)

ze_barplot <- barplot(NA, beside = TRUE, legend.text = T, col = alpha(cols, 0.2), ylim = c(0, lim), ylab = expression(Growth~(Delta*gr%*%day^-1)), xlab = "Environment", las = 1)
grid()

par(new = TRUE)
ze_barplot <- barplot(mean.gr, beside = TRUE, legend.text = T, col = alpha(cols, 0.5) , ylim = c(0, lim), ylab = expression(Growth~(Delta*gr%*%day^-1)), xlab = "Environment", las = 1)
error.bar(ze_barplot, mean.gr, stdev, col = "black")
box()
mtext('B', at = -0.45, line = -0.5)

##dev.off()


Violin plots

##png("figure2-revised.png", height = 7, width = 7, units = "in", res = 360)

layout(matrix(c(0, 1, 1, 0,
                0, 1, 1, 0,
                0, 2, 2, 0,
                0, 2, 2, 0), nrow = 4, ncol = 4, byrow = TRUE))

## Basic boxplot

par(mgp = c(2, 1, 0))

vioplot(prop_area_visited ~ strain + env, data = vid.dat, border = NA, method = "jitter", side = "right", ylab = expression(paste("Proportion of area visited")~(mm^2)), xlab = "Environment", col = "white", las = 1, xaxt = 'n', ylim = c(0, 0.45))

grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)

par(new = TRUE)

vioplot(prop_area_visited ~ strain + env, data = vid.dat, border = NA, method = "jitter", side = "right", ylab = "", xlab = "", col = c(alpha("purple", 0.2), alpha("orange", 0.2)), las = 1, xaxt = 'n', ylim = c(0, 0.45))
axis(1, at = c(1.5, 3.5), labels = c('Uniform', 'Patchy'), cex = 0.8)

legend('topleft', legend = c('Rover', 'Sitter'), cex = 0.8, fill = c(alpha("purple", 0.2), alpha("orange", 0.2)), bty = 'n')

stripchart(prop_area_visited ~ strain + env, data = vid.dat, vertical = TRUE, method = "jitter", add = TRUE, pch = 20, col = c(alpha("purple", 0.3), alpha("orange", 0.5)), cex = 1.3)

mtext("A", side = 2, at = 0.5, line = 3, las = 1, font = 1)


## Basic boxplot

par(mgp = c(3, 1, 0))

vioplot(growth ~ strain + env, data = data, border = NA, method = "jitter", side = "right", ylab = expression(Growth~(Delta*gr%*%day^-1)), xlab = "Environment", col = "white", las = 1, xaxt = 'n')

grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)

par(new = TRUE)

vioplot(growth ~ strain + env, data = data, border = NA, method = "jitter", side = "right", ylab = "", xlab = "", col = c(alpha("purple", 0.2), alpha("orange", 0.2)), las = 1, xaxt = 'n')
axis(1, at = c(1.5, 3.5), labels = c('Uniform', 'Patchy'), cex = 0.8)

stripchart(growth ~ strain + env, data = data, vertical = TRUE, method = "jitter", add = TRUE, pch = 20, col = c(alpha("purple", 0.3), alpha("orange", 0.5)), cex = 1.3)

mtext("B", side = 2, at = 0.00096, line = 3, las = 1, font = 1)

legend('topright', legend = c('Rover', 'Sitter'), cex = 0.8, fill = c(alpha("purple", 0.2), alpha("orange", 0.2)), bty = 'n')

##dev.off()